home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
FM Towns: Free Software Collection 4
/
FM Towns Free Software Collection 4 - Disc 1.iso
/
t_os
/
pi_v2
/
pi.c
next >
Wrap
Text File
|
1991-10-18
|
33KB
|
946 lines
/* ANSI-C(?)によるπ計算 VERSION 2.0 */
/* Copyright(c) Daisuke Takahashi 1991 */
/* 任意の瞬間にキー入力によって計算を中断したいのですが */
/* High-Cで signal.h が使えないので、一定時間毎にセーブさせています */
/* その代わり、DOS.H未使用の為、殆どのCでコンパイルできます(^_^) */
#include <stdio.h>
#include <stdlib.h>
#include <ctype.h>
#include <time.h>
#include <string.h>
/* ----------- definitions ------------ */
#define TEN8 100000000L /* 1億 10^8 */
#define MAX 4190000000L /* max of USL - TEN8 */
#define ON (1)
#define OFF (0)
#define SET (2)
#define LAP (1)
#define DONE (0)
#define DATA (1)
#define PLUS (1)
#define MINUS (-1)
#define NUL (\0)
#define USL unsigned long /* typedef にすべきでしょうが(^_^;) */
#define RATE (int)( 1000 * (long)ex.head / (long)ex.max )
/* ---------------- prototypes ------------------ */
void set_ex(int,char**); /* 引数→変数 */
void get_memory(void);
void back_memory(void); /* ワーク返還 */
void formal(void); /* 計算公式 */
void arctan(int,int,int); /* arctan設定 */
long bench_mark(long); /* 速度測定 */
void message(int,int,int); /* 計算量表示 */
void c_std(int,USL,USL,long*,long*); /* 41以下 */
void c_low(int,USL,USL,long*,long*); /* 6.4万以下 */
void c_mid(int,USL,USL,long*,long*); /* 200万以下 */
void c_high(int,USL,USL,long*,long*); /* 256万以下 */
void c_huge(int,USL,USL,long*,long*); /* 1600万以下 */
void regular(int); /* 正規化 */
void load(char *); /* 経過ロード */
void save(void); /* 経過セーブ */
void show_ans(void); /* 結果の表示 */
void time_set(void); /* 開始時刻記録 */
void time_end(void); /* 計算時間更新 */
void usage(int); /* 使い方 */
/* ------------- structs --------- */
struct {
/* ファイルセーブする変数 */
int ver; /* Version */
long sec; /* 計算時間 */
long cent; /* 1/100秒 */
long keta; /* 計算桁数 */
int max; /* 配列数 */
int type; /* 計算公式 */
int head; /* 動的先頭 */
int tail; /* 動的最終 */
USL save; /* 退避flag */
USL con; /* */
USL var; /* */
/* セーブ必要なし */
long *arc; /* 計算作業領域 */
long *ans; /* 計算結果領域 */
int load; /* file読込済 */
time_t s_sec; /* 計算開始時刻 */
clock_t s_cent; /* 〃1/100 秒 */
} ex;
int main( int argc, char **argv ) {
set_ex( argc, argv ); /* 引き数から変数設定 */
formal(); /* 計算 */
show_ans(); /* 結果と時間を表示 */
back_memory(); /* ワーク返還 */
return 0; /* 正常終了 */
}
void set_ex( int argc, char **argv ) {
char *buf;
if (argc == 1) /* オプション無しなら */
usage( !DATA ); /* 普通説明表示 */
if (strchr( argv[1], '.' ) != NULL) /* ファイル名であれば */
load( argv[1] ); /* 途中結果を読み込む */
else if ( !(isdigit( *argv[1] ) ) ) /* 無効オプションなら */
usage( DATA ); /* データ説明 */
else { /* 桁数指定オプション */
ex.ver = 0; /* On memory Version */
ex.sec = 0; /* 計算時間(秒単位) */
ex.cent = 0; /* 計算時間(0.01秒単位 */
ex.keta = atol( argv[1] ); /* 計算桁を代入 */
ex.max = (int)(ex.keta / 8 + 2); /* 配列数 */
ex.save = 0; /* セーブモードでない */
get_memory(); /* 必要メモリを確保 */
}
if (argc == 3) {
buf = argv[2];
while ( *buf ) {
if ( isdigit(*buf) ) /* 数字の場合は */
ex.save = strtol( buf, &buf, 10 ); /* 分単位のセーブ時間 */
else {
if (strchr( "gmrsGMRS", *buf ) != NULL && ex.load == OFF)
ex.type = *buf; /* ロード後の公式変更は不可 */
buf++; /* 計算公式を代入 */
}
}
}
}
void get_memory( void ) { /* ワークを確保 */
ex.arc = (long *)calloc( (size_t)ex.max, (size_t)sizeof(long) );
ex.ans = (long *)calloc( (size_t)ex.max, (size_t)sizeof(long) );
if (ex.ans == NULL) {
fputs( "メモリーが足らんぜよ(;_;) \n", stderr );
exit( 1 );
}
}
void back_memory( void ) { /* ワークを返還 */
free( ex.arc );
free( ex.ans );
}
void formal( void ) { /* 計算公式の値をセット */
switch ( tolower( ex.type ) ) {
case 'g': /* ガウスの式 */
/* 48*arctan(1/18) + 32*arctan(1/57) - 20*arctan(1/239) */
arctan( PLUS, 48, 18 );
arctan( PLUS, 32, 57 );
arctan( MINUS, 20, 239 );
break;
case 'r': /* ラザフォードの式 */
/* 16*arctan(1/5) - 4*arctan(1/70) + 4*arctan(1/99) */
arctan( PLUS, 16, 5 );
arctan( MINUS, 4, 70 );
arctan( PLUS, 4, 99 );
break;
case 's': /* シュテルマーの式 */
/* 24*arctan(1/8) + 8*arctan(1/57) + 4*arctan(1/239) */
arctan( PLUS, 24, 8 );
arctan( PLUS, 8, 57 );
arctan( PLUS, 4, 239 );
break;
default: /* デフォルトは最速マチン式 */
case 'm': /* マチンの式 */
/* 16*arctan(1/5) - 4*arctan(1/239) */
arctan( PLUS, 16, 5 );
arctan( MINUS, 4, 239 );
break;
}
if (ex.load == ON)
fputs( "Not calculated !!!\n", stderr );
}
void arctan( int sign, int mag, int base ) {
/* arctan(1/base) = 1/base - 1/(3*base^3) + 1/(5*base^5) - */
time_t cur_time, pre_time;
long min_one; /* 1分相当の var */
message( SET, sign, SET ); /* 3番目はダミー:計算式の符号 */
if (ex.load == OFF ) { /* 新規計算開始 */
*ex.arc = (long)base * (long)mag; /* 初期値設定 */
ex.con = (long)base * (long)base;
ex.head = 0;
ex.var = 1;
if ( base == 5 || base == 8 ) /* 10進数で有限小数なら */
ex.tail = 1; /* 最終位置は動的に変化 */
else /* 10進数で循環小数なら最終位置は物理的最終 */
ex.tail = ex.max;
} else { /* 計算再開 */
if (ex.con != (long)base * (long)base) {
message( DONE, mag, base );
return; /* この項は計算済 */
} else
ex.load = OFF; /* この項から計算再開:ロードフラグクリア */
}
message( LAP, mag, base );
time_set(); /* ここから計算時間測定開始 */
time( &pre_time ); /* 表示間隔初期値 */
min_one = bench_mark( 0 ); /* 時間測定初期化 */
while ( (ex.head < ex.max) || *(ex.arc + ex.max) ) {
USL check; /* どのcalc_*_? を使うかの判定 */
if ( ex.var > 64000L && (TEN8 % ex.var) < (MAX / ex.var) ) {
check = 60000L; /* enough under 64000L */
/* ans_r * T8_v_r < 41億 なら calc_low() で計算 */
} else
check = ex.var;
if (check < 42)
c_std( sign, ex.con, ex.var, ex.arc, ex.ans );
else if (check < 64000L)
c_low( sign, ex.con, ex.var, ex.arc, ex.ans );
else if (check < 2030000L)
c_mid( sign, ex.con, ex.var, ex.arc, ex.ans );
else if (check < 2560000L)
c_high( sign, ex.con, ex.var, ex.arc, ex.ans );
else /* check < 1600万 */
c_huge( sign, ex.con, ex.var, ex.arc, ex.ans );
sign *= MINUS; /* 符号逆転 */
if ( ex.tail < ex.max && *(ex.arc + ex.tail) )
ex.tail++;
if ( !*(ex.arc + ex.head) && ex.head < ex.max )
ex.head++; /* 有効数字の先頭に移動 */
ex.var += 2; /* 次の展開項を計算 */
if ( ex.var % 80 == 1 ) { /* ex.var 80 毎に */
regular( OFF ); /* ans[]を正規表現化(高速モード) */
if (min_one < 80) /* 1分相当値未計算なら */
min_one = bench_mark( min_one ); /* 計算する */
if ( (min_one >= 80) && ex.var % min_one == 1 ) { /* 1分経過なら */
time( &cur_time ); /* 現在時刻の取得 */
if (ex.save) { /* セーブ指定有りで時間なら */
if ( (long)difftime( cur_time, ex.s_sec ) > 60 * ex.save )
save(); /* セーブするだよん♪(^0^) */
}
if ( (long)difftime( cur_time, pre_time ) <= 40 ) {
min_one *= 2; /* 表示間隔が短くなったら倍にする */
}
pre_time = cur_time;
message( LAP, mag, base ); /* 表示時 Ctrl-C有効 */
}
}
}
regular( ON );
time_end(); /* 次のarctan()の為 */
message( DONE, mag, base );
if (ex.save) /* 一つのarctan(1/n)終了時も */
save(); /* saveする */
}
long bench_mark( long min_one ) { /* 1 分相当の ex.var 差を返す */
time_t cur_time;
long sec;
static USL st_var;
if (min_one == 0) { /* 初回計測? */
st_var = ex.var; /* 計測開始時刻 */
return 1; /* 初回時は計測 */
}
if (min_one >= 2) /* まだまだ時間がある */
return min_one - 1;
time( &cur_time );
sec = (long)difftime( cur_time, ex.s_sec );
if (sec >= 60)
return ex.var - st_var; /* 1分相当値:多桁のロスを防ぐ為80加算 */
else if (sec >= 30) /* 80加算の理由は上と同じ */
return (ex.var - st_var) *60 /sec /80 *80; /* 1分近似値 */
else /* 30秒以上余裕がある場合は */
return (60 - sec) / (sec + 1); /* 残り回数の推算 */
}
void message( int mode, int mag, int base ) {
static int sign;
if (mode == SET)
sign = mag; /* 計算式全体のフラグ */
if (mode == LAP)
fprintf( stderr, "%darctan(1/%d) %2d.%d%%\r",
sign * mag, base, RATE/10, RATE%10 );
if (mode == DONE)
fprintf( stderr, "%darctan(1/%d) completed !!!\n", sign * mag, base );
}
void c_std( int sign, USL con, USL var, long *arc, long *ans ) {
/* var must be under 42 */
int i, j;
USL co, va, buf, arc_r, ans_r;
USL con_q = TEN8 / con; /* arctan(1/5)以外 */
USL con_r = TEN8 % con; /* arctan(1/5)以外 */
if (sign == PLUS) {
if (con <= 40) { /* arctan(1/5)のみ */
for (co = con, va = var, arc_r = 0, ans_r = 0,
i = ex.head, j = ex.tail; i <= j; i++) {
buf = *(arc+i) + arc_r * TEN8;
*(arc+i) = buf / co;
arc_r = buf % co;
buf = *(arc+i) + ans_r * TEN8;
*(ans+i) += buf / va;
ans_r = buf % va;
}
for (j = ex.max; i <= j; i++) {
buf = ans_r * TEN8;
*(ans+i) += buf / va;
ans_r = buf % va;
}
} else { /* arctan(1/5)以外 */
for (co = con, va = var, arc_r = 0, ans_r = 0,
i = ex.head, j = ex.tail; i <= j; i++) {
buf = *(arc+i) + arc_r * con_r;
*(arc+i) = buf / co + arc_r * con_q;
arc_r = buf % co;
buf = *(arc+i) + ans_r * TEN8;
*(ans+i) += buf / va;
ans_r = buf % va;
}
for (j = ex.max; i <= j; i++) {
buf = ans_r * TEN8;
*(ans+i) += buf / va;
ans_r = buf % va;
}
}
} else { /* sign == MINUS */
if (con <= 40) { /* arctan(1/5)のみ */
for (co = con, va = var, arc_r = 0, ans_r = 0,
i = ex.head, j = ex.tail; i <= j; i++) {
buf = *(arc+i) + arc_r * TEN8;
*(arc+i) = buf / co;
arc_r = buf % co;
buf = *(arc+i) + ans_r * TEN8;
*(ans+i) -= buf / va;
ans_r = buf % va;
}
for (j = ex.max; i <= j; i++) {
buf = ans_r * TEN8;
*(ans+i) -= buf / va;
ans_r = buf % va;
}
} else { /* arctan(1/5)以外 */
for (co = con, va = var, arc_r = 0, ans_r = 0,
i = ex.head, j = ex.tail; i <= j; i++) {
buf = *(arc+i) + arc_r * con_r;
*(arc+i) = buf / co + arc_r * con_q;
arc_r = buf % co;
buf = *(arc+i) + ans_r * TEN8;
*(ans+i) -= buf / va;
ans_r = buf % va;
}
for (j = ex.max; i <= j; i++) {
buf = ans_r * TEN8;
*(ans+i) -= buf / va;
ans_r = buf % va;
}
}
}
}
void c_low( int sign, USL con, USL var, long *arc, long *ans ) {
/* var must be under 64_000 */
int i, j;
USL co, va, buf, arc_r, ans_r;
USL con_q = TEN8 / con; /* arctan(1/5)以外 */
USL con_r = TEN8 % con; /* arctan(1/5)以外 */
USL var_q = TEN8 / var;
USL var_r = TEN8 % var;
if (sign == PLUS) {
if (con <= 40) { /* arctan(1/5)のみ */
for ( co = con, va = var, arc_r = 0, ans_r = 0,
i = ex.head, j = ex.tail; i <= j; i++) {
buf = *(arc+i) + arc_r * TEN8;
*(arc+i) = buf / co;
arc_r = buf % co;
buf = *(arc+i) + ans_r * var_r;
*(ans+i) += buf / va + ans_r * var_q;
ans_r = buf % va;
}
for (j = ex.max; i <= j; i++ ) {
buf = ans_r * var_r;
*(ans+i) += buf / va + ans_r * var_q;
ans_r = buf % va;
}
} else { /* arctan(1/5)以外 */
for (co = con, va = var, arc_r = 0, ans_r = 0,
i = ex.head, j = ex.tail; i <= j; i++) {
buf = *(arc+i) + arc_r * con_r;
*(arc+i) = buf / co + arc_r * con_q;
arc_r = buf % co;
buf = *(arc+i) + ans_r * var_r;
*(ans+i) += buf / va + ans_r * var_q;
ans_r = buf % va;
}
for (j = ex.max; i <= j; i++ ) {
buf = ans_r * var_r;
*(ans+i) += buf / va + ans_r * var_q;
ans_r = buf % va;
}
}
} else { /* sign == MINUS */
if (con <= 40) { /* arctan(1/5)のみ */
for (co = con, va = var, arc_r = 0, ans_r = 0,
i = ex.head, j = ex.tail; i <= j; i++) {
buf = *(arc+i) + arc_r * TEN8;
*(arc+i) = buf / co;
arc_r = buf % co;
buf = *(arc+i) + ans_r * var_r;
*(ans+i) -= buf / va + ans_r * var_q;
ans_r = buf % va;
}
for (j = ex.max; i <= j; i++ ) {
buf = ans_r * var_r;
*(ans+i) -= buf / va + ans_r * var_q;
ans_r = buf % va;
}
} else { /* arctan(1/5)以外 */
for (co = con, va = var, arc_r = 0, ans_r = 0,
i = ex.head, j = ex.tail; i <= j; i++) {
buf = *(arc+i) + arc_r * con_r;
*(arc+i) = buf / co + arc_r * con_q;
arc_r = buf % co;
buf = *(arc+i) + ans_r * var_r;
*(ans+i) -= buf / va + ans_r * var_q;
ans_r = buf % va;
}
for (j = ex.max; i <= j; i++ ) {
buf = ans_r * var_r;
*(ans+i) -= buf / va + ans_r * var_q;
ans_r = buf % va;
}
}
}
}
void c_mid( int sign, USL con, USL var, long *arc, long *ans ) {
/* var must be under 2_030_000 */
int i, j;
USL co, va, buf, buf2, quot;
USL arc_r, ans_r;
USL con_q = TEN8 / con; /* arctan(1/5)以外 */
USL con_r = TEN8 % con; /* arctan(1/5)以外 */
USL var_q = TEN8 / var;
USL var_r1 = (TEN8 % var) / 1024;
USL var_r2 = (TEN8 % var) % 1024;
if (sign == PLUS) {
if (con <= 40) { /* arctan(1/5)のみ */
for (co = con, va = var, arc_r = 0, ans_r = 0,
i = ex.head, j = ex.tail; i <= j; i++) {
buf = *(arc+i) + arc_r * TEN8;
arc_r = buf % co;
*(arc+i) = buf / co;
buf = ans_r * var_r1;
quot = ( (buf / va) << 10 ) + ans_r * var_q;
buf2 = ( (buf % va) << 10 ) + ans_r * var_r2 + *(arc+i);
*(ans+i) += buf2 / va + quot;
ans_r = buf2 % va;
}
for (j = ex.max; i <= j; i++) {
buf = ans_r * var_r1;
quot = ( (buf / va) << 10 ) + ans_r * var_q;
buf2 = ( (buf % va) << 10 ) + ans_r * var_r2;
*(ans+i) += buf2 / va + quot;
ans_r = buf2 % va;
}
} else { /* arctan(1/5)以外 */
for (co = con, va = var, arc_r = 0, ans_r = 0,
i = ex.head, j = ex.tail; i <= j; i++) {
buf = *(arc+i) + arc_r * con_r;
*(arc+i) = buf / co + arc_r * con_q;
arc_r = buf % co;
buf = ans_r * var_r1;
quot = ( (buf / va) << 10 ) + ans_r * var_q;
buf2 = ( (buf % va) << 10 ) + ans_r * var_r2 + *(arc+i);
*(ans+i) += buf2 / va + quot;
ans_r = buf2 % va;
}
for (j = ex.max; i <= j; i++) {
buf = ans_r * var_r1;
quot = ( (buf / va) << 10 ) + ans_r * var_q;
buf2 = ( (buf % va) << 10 ) + ans_r * var_r2;
*(ans+i) += buf2 / va + quot;
ans_r = buf2 % va;
}
}
} else { /* sign == MINUS */
if (con <= 40) { /* arctan(1/5)のみ */
for (co = con, va = var, arc_r = 0, ans_r = 0,
i = ex.head, j = ex.tail; i <= j; i++) {
buf = *(arc+i) + arc_r * TEN8;
arc_r = buf % co;
*(arc+i) = buf / co;
buf = ans_r * var_r1;
quot = ( (buf / va) << 10 ) + ans_r * var_q;
buf2 = ( (buf % va) << 10 ) + ans_r * var_r2 + *(arc+i);
*(ans+i) -= buf2 / va + quot;
ans_r = buf2 % va;
}
for (j = ex.max; i <= j; i++) {
buf = ans_r * var_r1;
quot = ( (buf / va) << 10 ) + ans_r * var_q;
buf2 = ( (buf % va) << 10 ) + ans_r * var_r2;
*(ans+i) -= buf2 / va + quot;
ans_r = buf2 % va;
}
} else { /* arctan(1/5)以外 */
for (co = con, va = var, arc_r = 0, ans_r = 0,
i = ex.head, j = ex.tail; i <= j; i++) {
buf = *(arc+i) + arc_r * con_r;
*(arc+i) = buf / co + arc_r * con_q;
arc_r = buf % co;
buf = ans_r * var_r1;
quot = ( (buf / va) << 10 ) + ans_r * var_q;
buf2 = ( (buf % va) << 10 ) + ans_r * var_r2 + *(arc+i);
*(ans+i) -= buf2 / va + quot;
ans_r = buf2 % va;
}
for (j = ex.max; i <= j; i++) {
buf = ans_r * var_r1;
quot = ( (buf / va) << 10 ) + ans_r * var_q;
buf2 = ( (buf % va) << 10 ) + ans_r * var_r2;
*(ans+i) -= buf2 / va + quot;
ans_r = buf2 % va;
}
}
}
}
void c_high( int sign, USL con, USL var, long *arc, long *ans ) {
/* var must be under 2_560_000 */
int i, j;
USL co, va, buf, buf2, quot;
USL arc_r, ans_r;
USL con_q = TEN8 / con; /* arctan(1/5)以外 */
USL con_r = TEN8 % con; /* arctan(1/5)以外 */
USL var_q = TEN8 / var;
USL var_r1 = (TEN8 % var) / 1600;
USL var_r2 = (TEN8 % var) % 1600;
if (sign == PLUS) {
if (con <= 40) {
for (co = con, va = var, arc_r = 0, ans_r = 0,
i = ex.head, j = ex.tail; i <= j; i++) {
buf = *(arc+i) + arc_r * TEN8;
*(arc+i) = buf / co;
arc_r = buf % co;
buf = ans_r * var_r1;
quot = ans_r * var_q;
quot += 1600 * (buf / va);
buf2 = 1600 * (buf % va) + ans_r * var_r2 + *(arc+i);
*(ans+i) += buf2 / va + quot;
ans_r = buf2 % va;
}
for (j = ex.max; i <= j; i++) {
buf = ans_r * var_r1;
quot = ans_r * var_q;
quot += 1600 * (buf / va);
buf2 = 1600 * (buf % va) + ans_r * var_r2;
*(ans+i) += buf2 / va + quot;
ans_r = buf2 % va;
}
} else {
for (co = con, va = var, arc_r = 0, ans_r = 0,
i = ex.head, j = ex.tail; i <= j; i++) {
buf = *(arc+i) + arc_r * con_r;
*(arc+i) = buf / co + arc_r * con_q;
arc_r = buf % co;
buf = ans_r * var_r1;
quot = ans_r * var_q;
quot += 1600*(buf / va);
buf2 = 1600*(buf % va) + ans_r * var_r2 + *(arc+i);
*(ans+i) += buf2 / va + quot;
ans_r = buf2 % va;
}
for (j = ex.max; i <= j; i++) {
buf = ans_r * var_r1;
quot = ans_r * var_q;
quot += 1600 * (buf / va);
buf2 = 1600 * (buf % va) + ans_r * var_r2;
*(ans+i) += buf2 / va + quot;
ans_r = buf2 % va;
}
}
} else {
if (con <= 40) {
for (co = con, va = var, arc_r = 0, ans_r = 0,
i = ex.head, j = ex.tail; i <= j; i++) {
buf = *(arc+i) + arc_r * TEN8;
*(arc+i) = buf / co;
arc_r = buf % co;
buf = ans_r * var_r1;
quot = ans_r * var_q;
quot += 1600 * (buf / va);
buf2 = 1600 * (buf % va) + ans_r * var_r2 + *(arc+i);
*(ans+i) -= buf2 / va + quot;
ans_r = buf2 % va;
}
for (j = ex.max; i <= j; i++) {
buf = ans_r * var_r1;
quot = ans_r * var_q;
quot += 1600 * (buf / va);
buf2 = 1600 * (buf % va) + ans_r * var_r2;
*(ans+i) -= buf2 / va + quot;
ans_r = buf2 % va;
}
} else {
for (co = con, va = var, arc_r = 0, ans_r = 0,
i = ex.head, j = ex.tail; i <= j; i++) {
buf = *(arc+i) + arc_r * con_r;
*(arc+i) = buf / co + arc_r * con_q;
arc_r = buf % co;
buf = ans_r * var_r1;
quot = ans_r * var_q;
quot += 1600 * (buf / va);
buf2 = 1600 * (buf % va) + ans_r * var_r2 + *(arc+ i);
*(ans+i) -= buf2 / va + quot;
ans_r = buf2 % va;
}
for (j = ex.max; i <= j; i++) {
buf = ans_r * var_r1;
quot = ans_r * var_q;
quot += 1600 * (buf / va);
buf2 = 1600 * (buf % va) + ans_r * var_r2;
*(ans+i) -= buf2 / va + quot;
ans_r = buf2 % va;
}
}
}
}
void c_huge( int sign, USL con, USL var, long *arc, long *ans ) {
/* var must be under 15_625_000 (=250^3) */
int i, j;
USL co, va, buf, buf2, quot;
USL arc_r, ans_r;
USL con_q = TEN8 / con; /* arctan(1/5)以外 */
USL con_r = TEN8 % con; /* arctan(1/5)以外 */
USL var_q = TEN8 / var;
USL var_r1 = (TEN8 % var) / 62500L;
USL var_r2 = ( (TEN8 % var) % 62500L ) / 250;
USL var_r3 = (TEN8 % var) % 250;
if (sign == PLUS) {
if (con <= 40) { /* arctan(1/5)のみ */
for (co = con, va = var, arc_r = 0, ans_r = 0,
i = ex.head, j = ex.tail; i <= j; i++) {
buf = *(arc+i) + arc_r * TEN8;
*(arc+i) = buf / co;
arc_r = buf % co;
buf = ans_r * var_r1;
quot = ans_r * var_q;
quot += 62500L * (buf / va);
buf2 = 250 * (buf % va) + ans_r * var_r2;
quot += 250 * (buf2 / va);
buf = 250 * (buf2 % va) + ans_r * var_r3 + *(arc+i);
*(ans+i) += buf / va + quot;
ans_r = buf % va;
}
for (j = ex.max; i <= j; i++) {
buf = ans_r * var_r1;
quot = ans_r * var_q;
quot += 62500L * (buf / va);
buf2 = 250 * (buf % va) + ans_r * var_r2;
quot += 250 * (buf2 / va);
buf = 250 * (buf2 % va) + ans_r * var_r3;
*(ans+i) += buf / va + quot;
ans_r = buf % va;
}
} else { /* arctan(1/5)以外 */
for (co = con, va = var, arc_r = 0, ans_r = 0,
i = ex.head, j = ex.tail; i <= j; i++) {
buf = *(arc+i) + arc_r * con_r;
*(arc+i) = buf / co + arc_r * con_q;
arc_r = buf % co;
buf = ans_r * var_r1;
quot = ans_r * var_q;
quot += 62500L * (buf / va);
buf2 = 250 * (buf % va) + ans_r * var_r2;
quot += 250 * (buf2 / va);
buf = 250 * (buf2 % va) + ans_r * var_r3 + *(arc+i);
*(ans+i) += buf / va + quot;
ans_r = buf % va;
}
for (j = ex.max; i <= j; i++) {
buf = ans_r * var_r1;
quot = ans_r * var_q;
quot += 62500L * (buf / va);
buf2 = 250 * (buf % va) + ans_r * var_r2;
quot += 250 * (buf2 / va);
buf = 250 * (buf2 % va) + ans_r * var_r3;
*(ans+i) += buf / va + quot;
ans_r = buf % va;
}
}
} else { /* sign == MINUS */
if (con <= 40) { /* arctan(1/5)のみ */
for (co = con, va = var, arc_r = 0, ans_r = 0,
i = ex.head, j = ex.tail; i <= j; i++) {
buf = *(arc+i) + arc_r * TEN8;
*(arc+i) = buf / co;
arc_r = buf % co;
buf = ans_r * var_r1;
quot = ans_r * var_q;
quot += 62500L * (buf / va);
buf2 = 250 * (buf % va) + ans_r * var_r2;
quot += 250 * (buf2 / va);
buf = 250 * (buf2 % va) + ans_r * var_r3 + *(arc+i);
*(ans+i) -= buf / va + quot;
ans_r = buf % va;
}
for (j = ex.max; i <= j; i++) {
buf = ans_r * var_r1;
quot = ans_r * var_q;
quot += 62500L * (buf / va);
buf2 = 250 * (buf % va) + ans_r * var_r2;
quot += 250 * (buf2 / va);
buf = 250 * (buf2 % va) + ans_r * var_r3;
*(ans+i) -= buf / va + quot;
ans_r = buf % va;
}
} else { /* arctan(1/5)以外 */
for (co = con, va = var, arc_r = 0, ans_r = 0,
i = ex.head, j = ex.tail; i <= j; i++) {
buf = *(arc+i) + arc_r * con_r;
*(arc+i) = buf / co + arc_r * con_q;
arc_r = buf % co;
buf = ans_r * var_r1;
quot = ans_r * var_q;
quot += 62500L * (buf / va);
buf2 = 250 * (buf % va) + ans_r * var_r2;
quot += 250 * (buf2 / va);
buf = 250 * (buf2 % va) + ans_r * var_r3 + *(arc+i);
*(ans+i) -= buf / va + quot;
ans_r = buf % va;
}
for (j = ex.max; i <= j; i++) {
buf = ans_r * var_r1;
quot = ans_r * var_q;
quot += 62500L * (buf / va);
buf2 = 250 * (buf % va) + ans_r * var_r2;
quot += 250 * (buf2 / va);
buf = 250 * (buf2 % va) + ans_r * var_r3;
*(ans+i) -= buf / va + quot;
ans_r = buf % va;
}
}
}
}
void regular( int mode ) { /* ans[]を0~TEN8-1の範囲に収める */
int i;
long buf, carry;
long *ans = ex.ans;
buf = 0;
carry = 0;
if (mode == OFF) { /* 計算中は高速化の為、負の処理をしない */
for (i = ex.max; i >= 0; i--) {
buf = *(ans + i) + carry;
*(ans + i) = buf % TEN8;
carry = buf / TEN8;
}
} else { /* 負の処理も行う */
for (i = ex.max; i >= 0; i--) {
buf = *(ans + i) + carry;
*(ans + i) = buf % TEN8;
carry = buf / TEN8;
if ( *(ans + i) < 0) { /* 負の% 演算の処理 */
*(ans + i) += TEN8;
carry--;
}
}
}
}
void load( char *fname ) {
FILE *fp;
char buf[80];
int m;
if ( ( fp = fopen(fname, "r") ) == NULL ) {
fputs( "file not found", stderr );
usage( !DATA );
}
ex.load = ON; /* ロードフラグON */
fgets( buf, 79, fp ); /* まずコメント行をダミーリード */
fgets( buf, 79, fp ); /* 次に配列に変数群を読み込む */
if ( sscanf( buf, "%d%ld%ld%ld%d%d%d%d%ld%ld%ld",
&ex.ver, /* Version */
&ex.sec, /* 計算時間 */
&ex.cent, /* 〃1/100 */
&ex.keta, /* 計算桁数 */
&ex.max, /* 静的最終 */
&ex.type, /* 計算公式 */
&ex.head, /* 動的先頭 */
&ex.tail, /* 動的最終 */
&ex.save, /* 退避flag */
&ex.con, /* base^2 */
&ex.var ) /* arc変数 */
!= 11 /* 変数が足りない */
|| (ex.max != (int)ex.keta/8 + 2) /* 関係式がおかしい */
|| ! (ex.save && ex.con && ex.var) /* どれかが0 である */
) { /* どれでも成立すれば異常 */
fclose( fp );
fputs( "Not regular file", stderr );
exit( 3 );
}
get_memory();
for ( m = 0; m <= ex.max; m++ ) {
fgets( buf, 79, fp );
if ( sscanf( buf, "%ld%ld", ex.ans+m, ex.arc+m ) != 2 ) {
fclose( fp );
fputs( "Not regular file\n", stderr );
exit( 3 );
}
}
fclose( fp );
}
void save( void ) {
FILE *fp;
int m;
time_t cur_time; /* 現在時刻記録用 */
char fname[13]; /* 8+1+3+NUL = 13 */
char buf[80];
time_end();
sprintf( fname, "%ld.cal", ex.keta );
if ( ( fp = fopen(fname, "w") ) == NULL ) {
time( &cur_time );
fprintf( stderr, "file can't be opened at %s", ctime( &cur_time ) );
time_set(); /* セーブ失敗でも計算は続行 */
return;
}
regular( ON ); /* 既にregular()しているのでこれはズルでない */
fputs( "PI.C(Ver 2.0)(Normal mode) ここはコメント行(79文字以下)\n", fp );
sprintf( buf, "%d %ld %ld %ld %d %d %d %d %ld %ld %ld\n",
ex.ver, /* Version */
ex.sec, /* 計算時間 */
ex.cent, /* 〃1/100 */
ex.keta, /* 計算桁数 */
ex.max, /* 静的最終 */
ex.type, /* 計算公式 */
ex.head, /* 動的先頭 */
ex.tail, /* 動的最終 */
ex.save, /* 退避flag */
ex.con, /* base^2 */
ex.var ); /* arc変数 */
fputs( buf, fp );
for ( m = 0; m <= ex.max; m++ ) {
sprintf( buf, "%ld %ld\n", *(ex.ans+m), *(ex.arc+m) );
fputs( buf, fp );
}
fclose( fp );
time_set(); /* 計算開始時刻を再セット */
}
void time_set(void) { /* 開始時間セット */
time( &ex.s_sec ); /* 計算開始時刻(秒) */
ex.s_cent = clock(); /* 〃1/100 秒 */
}
void time_end( void ) { /* 経過時間更新 */
long sec, cent;
time_t e_sec;
clock_t e_cent;
time( &e_sec );
e_cent = clock();
sec = (long)difftime( e_sec, ex.s_sec );
cent = ( e_cent - ex.s_cent + 86400 * (long)CLK_TCK )
% ( 86400 * (long)CLK_TCK );
sec /= 86400;
sec *= 86400; /* 一日未満の秒は一旦切捨て、 */
sec += (cent / (long)CLK_TCK); /* cent より求めて再加算する */
cent *= 100; /* 1/100 秒単位の数にする為 */
cent /= (long)CLK_TCK; /* これでなっているハズ */
cent = cent % 100; /* 1 秒未満 */
ex.sec += sec; /* 計算時間更新 */
ex.cent += cent; /* 〃 */
if (ex.cent >= 100) {
ex.cent -= 100;
ex.sec++;
}
ex.s_sec = e_sec; /* 再開時刻更新(time_setと同じ) */
ex.s_cent = e_cent; /* 〃1/100 秒 */
}
void show_ans( void ) {
int i = 0; /* 配列count */
int l = 0; /* 文字列位置 */
char buf[60];
char *buf_p = buf;
fprintf( stderr, "計算時間は%ld時間", ex.sec/3600 );
fprintf( stderr, "%0.2ld分%0.2ld.%0.2ld秒でした\n\n",
(ex.sec/60)%60, ex.sec%60, ex.cent );
printf("πの小数点以下%d(+α)桁の値です。\n", ex.keta );
printf( "00000001:%0.1ld.", *ex.ans );
for ( i = 1; i <= ex.max; i++) {
sprintf( buf_p, "%0.8ld\0", *(ex.ans + i) );
buf_p += 8;
if ( strlen(buf) >= 50 || i == ex.max ) {
for ( l = 0; l < 50 && buf[l]; l++) {
putchar( buf[l] );
if ( l % 10 == 9 )
putchar( ' ' );
}
if (l != 50)
break;
printf( "\n%0.8d: ", (i*8 / 10)*10 + 1 );
strcpy( buf, buf + 50 ); /* 50桁分詰める */
buf_p -= 50; /* 50桁分戻る */
}
}
printf( "\n計算時間は%ld時間", ex.sec/3600 );
printf( "%0.2ld分%0.2ld.%0.2ld秒でした\n",
(ex.sec/60)%60, ex.sec%60, ex.cent );
}
void usage( int mode ) {
if (mode == DATA) {
puts("TECNICAL DATA・・・・通常の使い方は A>RUN386 PI で表示\n");
puts("TOWNS 2F(初代) 1万桁 10万桁 100万桁");
puts("ウェイト無調整 96秒 165分 333時間");
puts("公式別時間 マチン ガウス シュテルマー ラザフォード");
puts(" 96秒 131秒 141秒 141秒");
puts("\nセーブファイルについて");
puts(" サイズ: 桁数×(1.5~2.1)倍");
puts(" ファイル名:セーブ時:桁数.CAL(.CALは自動的に付く)");
puts(" ロード時:拡張子も含め任意の名が可能");
puts("\n計算再開時にセーブ時間の再変更が可能");
puts("セーブ指定有りでは、各arctan()終了毎にもセーブする。");
puts("\n計算進行表示は約1分毎だが途中多少の変動あり");
puts("\nまた、第2オプションは複数記述が可能");
exit(0);
}
puts("RUN386 PI 100000 ←第1オプション:桁数 or ファイル名");
puts(" 100000.CAL 拡張子必要:ファイルを読み込んで計算再開");
puts(" (others) 無効オプション: Tecnical data 表示");
puts("RUN386 PI 100000 G20 ←第2オプション:公式 or セーブ時間間隔(分)");
puts(" g:ガウス m:マチン(デフォルト) r:ラザフォード s:シュテルマー");
puts(" 例はガウス公式で20分おきにセーブする設定");
puts("約1分毎の計算進行表示時、 Ctrl-C による中断が可能である");
exit(0);
}